pacman::p_load(sp, sf, spNetwork, tmap,tidyverse,spatstat,raster,maptools,kableExtra,plotly,ggthemes,onemapsgapi,devtools,rgdal)Airbnb_Project
Import packages
Data
Geospatial Data
Import Data
beijing_sf <- st_read("data/Geospatial/beijing.geojson") %>%
st_transform(crs=4555)Reading layer `beijing' from data source
`C:\Quanfang777\IS415-GAA\Project\data\Geospatial\beijing.geojson'
using driver `GeoJSON'
Simple feature collection with 16 features and 2 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 115.4172 ymin: 39.44247 xmax: 117.5071 ymax: 41.05861
Geodetic CRS: WGS 84Select Column for Analysis
beijing_sf <- select(beijing_sf, neighbourhood, geometry)
beijing_sfSimple feature collection with 16 features and 1 field
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 115.4172 ymin: 39.44247 xmax: 117.5071 ymax: 41.05861
Geodetic CRS: New Beijing
First 10 features:
neighbourhood geometry
1 东城区 MULTIPOLYGON (((116.4423 39...
2 西城区 MULTIPOLYGON (((116.3916 39...
3 昌平区 MULTIPOLYGON (((116.0427 40...
4 大兴区 / Daxing MULTIPOLYGON (((116.7347 39...
5 房山区 MULTIPOLYGON (((116.2466 39...
6 怀柔区 / Huairou MULTIPOLYGON (((116.279 40....
7 门头沟区 / Mentougou MULTIPOLYGON (((115.563 39....
8 密云县 / Miyun MULTIPOLYGON (((116.8826 40...
9 平谷区 / Pinggu MULTIPOLYGON (((117.3813 40...
10 延庆县 / Yanqing MULTIPOLYGON (((116.279 40....Check CRS Again
st_crs(beijing_sf)Coordinate Reference System:
User input: EPSG:4555
wkt:
GEOGCRS["New Beijing",
DATUM["New Beijing",
ELLIPSOID["Krassowsky 1940",6378245,298.3,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Topographic mapping."],
AREA["China - onshore."],
BBOX[18.11,73.62,53.56,134.77]],
ID["EPSG",4555]]Check if there is invalid geometry
length(which(st_is_valid(beijing_sf ) == FALSE))[1] 0Check Missing Value
beijing_sf[rowSums(is.na(beijing_sf))!=0,]Simple feature collection with 0 features and 1 field
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: New Beijing
[1] neighbourhood geometry
<0 rows> (or 0-length row.names)Aspatial Data
beforecovid <- read_csv("data/aspatial/beijing_beforecovid.csv")Warning: One or more parsing issues, call `problems()` on your data frame for details,
e.g.:
dat <- vroom(...)
problems(dat)Rows: 25921 Columns: 106
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (51): listing_url, last_scraped, name, summary, space, description, expe...
dbl (39): id, scrape_id, host_id, host_listings_count, host_total_listings_c...
lgl (16): thumbnail_url, medium_url, xl_picture_url, host_is_superhost, host...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.aftercovid <- read_csv("data/aspatial/beijing_aftercovid.csv")Rows: 6050 Columns: 18
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (5): name, host_name, neighbourhood, room_type, last_review
dbl (11): id, host_id, latitude, longitude, price, minimum_nights, number_of...
lgl (2): neighbourhood_group, license
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Select Column for Analysis
beforecovid <- subset(beforecovid, select=c(id, street,neighbourhood_cleansed,latitude, longitude))aftercovid <- subset(aftercovid, select=c(id,neighbourhood,latitude, longitude))Check Missing Value
sum(is.na(beforecovid$latitude))[1] 0sum(is.na(aftercovid$latitude))[1] 0sum(is.na(beforecovid$longitude))[1] 0sum(is.na(aftercovid$longitude))[1] 0Visualize the listings
beforecovid_sf <-
st_as_sf(beforecovid,
# our coordinates are the longitude and latitude
coords = c("longitude",
"latitude"),
# the geographical features are in longitude & latitude, in decimals
# as such, WGS84 is the most appropriate coordinates system
crs = 4326) %>%st_transform(crs=4555)
#afterwards, we transform it to SVY21, our desired CRSaftercovid_sf <- st_as_sf(aftercovid,
# our coordinates are the longitude and latitude
coords = c("longitude",
"latitude"),
# the geographical features are in longitude & latitude, in decimals
# as such, WGS84 is the most appropriate coordinates system
crs = 4326) %>% st_transform(crs=4555)The listing before covid
tmap_mode("plot")+
qtm(beijing_sf) +
tm_shape(beforecovid_sf)+
tm_dots()tmap mode set to plotting
We could notice that there are quite a lot of listings locate at the central area of Beijing
The listing after covid
tmap_mode("plot")+
qtm(beijing_sf) +
tm_shape(aftercovid_sf)+
tm_dots()tmap mode set to plotting
After covid, the listings has reduced significantly
Network Analysis of Road
Import the Road Data of Beijing and name it as ‘network’
network <- st_read(dsn="data/Geospatial/Network/shape",
layer="roads")%>%
st_transform(crs=4555)%>% filter(type %in% c("residential", "footway","path"))Reading layer `roads' from data source
`C:\Quanfang777\IS415-GAA\Project\data\Geospatial\Network\shape'
using driver `ESRI Shapefile'
Simple feature collection with 98252 features and 7 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 116.08 ymin: 39.68 xmax: 116.77 ymax: 40.17999
Geodetic CRS: WGS 84network_sf <- select(network,osm_id,type, geometry)
network_sfSimple feature collection with 40274 features and 2 fields
Geometry type: LINESTRING
Dimension: XY
Bounding box: xmin: 116.08 ymin: 39.68 xmax: 116.7699 ymax: 40.17999
Geodetic CRS: New Beijing
First 10 features:
osm_id type geometry
1 4822342 residential LINESTRING (116.3972 39.993...
2 4964849 residential LINESTRING (116.4645 39.953...
3 5168733 footway LINESTRING (116.4561 39.942...
4 5168864 residential LINESTRING (116.4559 39.945...
5 5212775 footway LINESTRING (116.4561 39.945...
6 5217979 footway LINESTRING (116.4561 39.937...
7 9701854 residential LINESTRING (116.4584 39.952...
8 9779360 residential LINESTRING (116.401 39.9116...
9 9932203 residential LINESTRING (116.4611 39.932...
10 14347213 residential LINESTRING (116.4424 39.935...Check the CRS
st_crs(network_sf)Coordinate Reference System:
User input: EPSG:4555
wkt:
GEOGCRS["New Beijing",
DATUM["New Beijing",
ELLIPSOID["Krassowsky 1940",6378245,298.3,
LENGTHUNIT["metre",1]]],
PRIMEM["Greenwich",0,
ANGLEUNIT["degree",0.0174532925199433]],
CS[ellipsoidal,2],
AXIS["geodetic latitude (Lat)",north,
ORDER[1],
ANGLEUNIT["degree",0.0174532925199433]],
AXIS["geodetic longitude (Lon)",east,
ORDER[2],
ANGLEUNIT["degree",0.0174532925199433]],
USAGE[
SCOPE["Topographic mapping."],
AREA["China - onshore."],
BBOX[18.11,73.62,53.56,134.77]],
ID["EPSG",4555]]Check invalid geometry and missing value
length(which(st_is_valid(network_sf) == FALSE))[1] 0network_sf[rowSums(is.na(network_sf))!=0,]Simple feature collection with 0 features and 2 fields
Bounding box: xmin: NA ymin: NA xmax: NA ymax: NA
Geodetic CRS: New Beijing
[1] osm_id type geometry
<0 rows> (or 0-length row.names)Focus on Three Study Area
1. Shijingshan District(石景山区)
shijingshan_sf <- beijing_sf %>% filter(neighbourhood == "石景山区")plot(shijingshan_sf)
Only select the road that within Shijingshan District
network_shijingshan <- subset(network_sf, lengths(st_intersects(network_sf, shijingshan_sf))!=0,)Only Select the listings that within Shijingshan District
shijingshan_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, shijingshan_sf))!=0,)
shijingshan_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, shijingshan_sf))!=0,)Visualize the Listings and the Network in Shijingshan Before Covid
tmap_mode("plot")tmap mode set to plottingtm_shape(shijingshan_sf) +
tm_polygons() +
tm_shape(network_shijingshan) +
tm_lines() +
tm_shape(shijingshan_beforecovidlisting) +
tm_dots(size = 0.01,
col = "blue",
border.col="black",
border.lwd=0.5)
Visualize the Listings and the Network in Shijingshan After Covid
tmap_mode("plot")tmap mode set to plottingtm_shape(shijingshan_sf) +
tm_polygons() +
tm_shape(network_shijingshan) +
tm_lines() +
tm_shape(shijingshan_aftercovidlisting) +
tm_dots(size = 0.01,
col = "blue",
border.col="black",
border.lwd=0.5)
Assign CRS
crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_shijingshan<- st_transform(network_shijingshan, crs)
shijingshan_beforecovidlisting<- st_transform(shijingshan_beforecovidlisting, crs)
shijingshan_aftercovidlisting<- st_transform(shijingshan_aftercovidlisting, crs)Network Constrained KDE (NetKDE) Analysis
Preparing the lixels objects
Before computing NetKDE, the SpatialLines object need to be cut into lixels with a specified minimal distance
lixels <- lixelize_lines(network_shijingshan,
700,
mindist = 350)Generating line centre points
Next, lines_center() of spNetwork will be used to generate a SpatialPointsDataFrame (i.e. samples) with line centre points as shown in the code chunk below.
samples <- lines_center(lixels)Performing NetKDE for Beijing Before Covid
densities <- nkde(network_shijingshan,
events = shijingshan_beforecovidlisting,
w = rep(1,nrow(shijingshan_beforecovidlisting)),
samples = samples,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 10,
sparse = TRUE,
verbose = FALSE)samples$density <- densities
lixels$density <- densitiessamples$density <- samples$density*1000
lixels$density <- lixels$density*1000tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels)+
tm_lines(col="density")+
tm_shape(shijingshan_beforecovidlisting)+
tm_dots()Network Constrained G- and K-Function Analysis
Performing NetKDE for Beijing After Covid
densities1 <- nkde(network_shijingshan,
events = shijingshan_aftercovidlisting,
w = rep(1,nrow(shijingshan_aftercovidlisting)),
samples = samples,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg =10,
sparse = TRUE,
verbose = FALSE)samples$density_after <- densities1
lixels$density_after<- densities1
samples$density_after <- samples$density_after*1000
lixels$density_after <- lixels$density_after*1000tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels)+
tm_lines(col="density_after")+
tm_shape(shijingshan_aftercovidlisting)+
tm_dots()2. Haidian 海淀
haidian_sf <- beijing_sf %>% filter(neighbourhood == "海淀区")%>%
st_transform(crs=4555)plot(haidian_sf)
Only select the road that within Haidian District
network_haidian <- subset(network_sf, lengths(st_intersects(network_sf, haidian_sf ))!=0,)%>%
st_transform(crs=4555)haidian_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, haidian_sf))!=0,) %>%
st_transform(crs=4555)
haidian_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, haidian_sf))!=0,)%>%
st_transform(crs=4555)Visualize the Listings and the Network in Haidian Before Covid
tmap_mode("plot")tmap mode set to plottingtm_shape(haidian_sf) +
tm_polygons() +
tm_shape(network_haidian) +
tm_lines() +
tm_shape(haidian_beforecovidlisting) +
tm_dots(size = 0.01,
col = "blue",
border.col="black",
border.lwd=0.5)
Visualize the Listings and the Network in Haidian After Covid
tmap_mode("plot")tmap mode set to plottingtm_shape(haidian_sf) +
tm_polygons() +
tm_shape(network_haidian) +
tm_lines() +
tm_shape(haidian_aftercovidlisting) +
tm_dots(size = 0.01,
col = "blue",
border.col="black",
border.lwd=0.5)
crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_haidian<- st_transform(network_haidian, crs)
haidian_beforecovidlisting<- st_transform(haidian_beforecovidlisting, crs)
haidian_aftercovidlisting<- st_transform(haidian_aftercovidlisting, crs)Network Constrained KDE (NetKDE) Analysis
lixels_haidian <- lixelize_lines(network_haidian,
700,
mindist = 350)samples_haidian <- lines_center(lixels_haidian)samples_haidian<- st_transform(samples_haidian, crs)densities_haidian <- nkde(network_haidian,
events = haidian_beforecovidlisting,
w = rep(1,nrow(haidian_beforecovidlisting)),
samples = samples_haidian,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 10,
sparse = TRUE,
verbose = FALSE)samples_haidian$density <- densities_haidian
lixels_haidian$density <- densities_haidiansamples_haidian$density <- samples_haidian$density*1000
lixels_haidian$density <- lixels_haidian$density*1000The listings before covid
tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels_haidian)+
tm_lines(col="density")+
tm_shape(haidian_beforecovidlisting)+
tm_dots()After Covid Visualization
densities_haidian <- nkde(network_haidian,
events = haidian_aftercovidlisting,
w = rep(1,nrow(haidian_aftercovidlisting)),
samples = samples_haidian,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 10,
sparse = TRUE,
verbose = FALSE)samples_haidian$density_after <- densities_haidian
lixels_haidian$density_after <- densities_haidian
samples_haidian$density_after <- samples_haidian$density_after*1000
lixels_haidian$density_after <- lixels_haidian$density_after*1000tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels_haidian)+
tm_lines(col="density")+
tm_shape(haidian_aftercovidlisting)+
tm_dots()FENGTAI 丰台
fengtai_sf <- beijing_sf %>% filter(neighbourhood == "丰台区 / Fengtai")%>%
st_transform(crs=4555)
plot(fengtai_sf)
network_fengtai <- subset(network_sf, lengths(st_intersects(network_sf, fengtai_sf ))!=0,)%>%
st_transform(crs=4555)
plot(network_fengtai)
Only select the road that within Fengtai District
fengtai_beforecovidlisting <- subset(beforecovid_sf, lengths(st_intersects(beforecovid_sf, fengtai_sf))!=0,) %>%
st_transform(crs=4555)
fengtai_aftercovidlisting <- subset(aftercovid_sf, lengths(st_intersects(aftercovid_sf, fengtai_sf))!=0,)%>%
st_transform(crs=4555)Visualize the Listings and the Network in Haidian Before Covid
tmap_mode("plot")tmap mode set to plottingtm_shape(fengtai_sf) +
tm_polygons() +
tm_shape(network_fengtai) +
tm_lines() +
tm_shape(fengtai_beforecovidlisting) +
tm_dots(size = 0.01,
col = "blue",
border.col="black",
border.lwd=0.5)
Visualize the Listings and the Network in Haidian Before Covid
tmap_mode("plot")tmap mode set to plottingtm_shape(fengtai_sf) +
tm_polygons() +
tm_shape(network_fengtai) +
tm_lines() +
tm_shape(fengtai_aftercovidlisting) +
tm_dots(size = 0.01,
col = "blue",
border.col="black",
border.lwd=0.5)
Network Constrained KDE (NetKDE) Analysis
Preparing the lixels objects
crs <- "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
network_fengtai<- st_transform(network_fengtai, crs)
fengtai_beforecovidlisting<- st_transform(fengtai_beforecovidlisting, crs)
fengtai_aftercovidlisting<- st_transform(fengtai_aftercovidlisting, crs)lixels_fengtai <- lixelize_lines(network_fengtai,
700,
mindist = 350)samples_fengtai <- lines_center(lixels_fengtai)densities_fengtai <- nkde(network_fengtai,
events = fengtai_beforecovidlisting,
w = rep(1,nrow(fengtai_beforecovidlisting)),
samples = samples_fengtai,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 10,
sparse = TRUE,
verbose = FALSE)samples_fengtai$density <- densities_fengtai
lixels_fengtai$density <- densities_fengtai
samples_fengtai$density <- samples_fengtai$density*1000
lixels_fengtai$density <- lixels_fengtai$density*1000tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels_fengtai)+
tm_lines(col="density")+
tm_shape(fengtai_beforecovidlisting)+
tm_dots()After covid
densities_fengtai_after <- nkde(network_fengtai,
events = fengtai_aftercovidlisting,
w = rep(1,nrow(fengtai_aftercovidlisting)),
samples = samples_fengtai,
kernel_name = "quartic",
bw = 300,
div= "bw",
method = "simple",
digits = 1,
tol = 1,
grid_shape = c(1,1),
max_depth = 8,
agg = 10,
sparse = TRUE,
verbose = FALSE)samples_fengtai$density_after <- densities_fengtai_after
lixels_fengtai$density_after <- densities_fengtai_after
samples_fengtai$density_after <- samples_fengtai$density_after*1000
lixels_fengtai$density_after <- lixels_fengtai$density_after*1000tmap_mode('view')tmap mode set to interactive viewingtm_shape(lixels_fengtai)+
tm_lines(col="density_after")+
tm_shape(fengtai_aftercovidlisting)+
tm_dots()The graphs reveal that road segments (darker orange color) are the part of network with relatively higher density of the listings. Those roads are located at/near the edge of the district The patterns are applicable to all of the three districts (Shijingshan, Fengtai, Haidian)
We could observe that the dense central area has been affected the most, the regions that remain unaffected are located at the edge of the district, showing that the lockdown measures has greatly impact the most of listings, especially those in the center of the Beijing where more tourists attraction are located, those areas are also the area with strict lockdown policy implementation
Insights Gained
The measures implemented to curb the spread of the pandemic had significant geographical effects. The closure or restricted access to public attractions, mostly located in central areas, greatly impacted the home-sharing and accommodation businesses like Airbnb. Risk management in the hospitality sector should be aware of changes in geographical spread of lodging facilities and consider how to strategise and optimise future locations of such facilities to better serve the needs of the market.